Krill Length 48.1
Mauricio Mardones
09 March, 2023
rm(list = ls())Background
The following document intends to carry out a complementary methodological analysis to correlate environmental variables with the population dynamics of krill (Euphausia superba).
Objective
Once the correlation and effects on the population and/or fishing indicators on krill have been verified, this analysis aims to have a time series of the environmental variable to incorporate into the stock assessment process. Similar work in Wang et al. (2021) but with a longest fishery history.
library(reshape2)
library(tidyverse)
library(plyr)
library(lubridate)
library(raster)
library(sf)
library(CCAMLRGIS)
library(here)
library(easystats) # analisis estadisticos post
library(ncdf4)
library(readr)
library(data.table)
library(ggridges)
library(patchwork)
library(knitr)
library(kableExtra)Length structure krill 48 subarea
Another important piece of information for a stock evaluation refers to the biological components such as average sizes and weights across areas and years. To do this, we will explore the biological data and prepare the output to pass it to SS3.
the idea is to correlate mean lengths with SIC.
Data exploratory analysis
The object ohbio2 come from data exploration analysis in
data request CCAMLR data. This objetc have bio information from
krill.
#datos entregados por secretaria de CCMLAR
metadata <- load("~/DOCAS/Data/565_C1_KRI_2021-10-01/DATA_PRODUCT_565.RData")
# Data procesada por MMardones
#ohbio <- load("DataLengthKrill.RData")
#ohbio
metadata## [1] "C1" "OBS_HAUL_BIOLOGY" "METADATA"#son lo mismo#cargo objeto
meta <- get("METADATA")
c1 <- get("C1")
ohbio <- get("OBS_HAUL_BIOLOGY")
names(ohbio)## [1] "c1_id" "obs_haul_id" "obs_logbook_id"
## [4] "haul_number" "taxon_code" "taxon_scientific_name"
## [7] "taxon_family" "maturity_stage" "sex_code"
## [10] "length_total_cm" "greenweight_kg"dim(c1)## [1] 360439 35dim(ohbio)## [1] 2107887 11Join data set with master as c1 set. This join is
trought obs_haul_id variable
ohbio2 <- left_join(c1, ohbio, by="obs_haul_id")
dim(ohbio2)## [1] 2443773 45Firsts glance. Test how many register have by year. In this case,
length_total_cm by season ccamlr.
same exercise in date period date_catchperiod_start
ohbio3 <- ohbio2 %>%
mutate(Year = year(date_catchperiod_start),
Month = month(date_catchperiod_start),
Day = day(date_catchperiod_start))names(ohbio3)## [1] "c1_id.x" "obs_haul_id"
## [3] "obs_logbook_id.x" "obs_haul_number"
## [5] "haul_number.x" "vessel_name"
## [7] "vessel_nationality_code" "fishing_purpose_code"
## [9] "season_ccamlr" "target_species"
## [11] "asd_code" "trawl_technique"
## [13] "catchperiod_code" "date_catchperiod_start"
## [15] "datetime_set_start" "datetime_set_end"
## [17] "datetime_haul_start" "datetime_haul_end"
## [19] "datetime_timezone" "depth_gear_set_end_m"
## [21] "depth_gear_haul_start_m" "depth_bottom_set_end_m"
## [23] "depth_bottom_haul_start_m" "latitude_set_end"
## [25] "longitude_set_end" "latitude_haul_start"
## [27] "longitude_haul_start" "gear_type_code"
## [29] "gear_type" "mesh_code"
## [31] "trawl_net_number" "notes"
## [33] "trawl_duration_depth_h" "trawl_duration_total_h"
## [35] "krill_greenweight_kg" "c1_id.y"
## [37] "obs_logbook_id.y" "haul_number.y"
## [39] "taxon_code" "taxon_scientific_name"
## [41] "taxon_family" "maturity_stage"
## [43] "sex_code" "length_total_cm"
## [45] "greenweight_kg" "Year"
## [47] "Month" "Day"ohbio4 <- ohbio3 %>%
dplyr::select(7, 9, 11, 12, 14, 24, 25, 29, 42, 44, 46, 47)
names(ohbio4)## [1] "vessel_nationality_code" "season_ccamlr"
## [3] "asd_code" "trawl_technique"
## [5] "date_catchperiod_start" "latitude_set_end"
## [7] "longitude_set_end" "gear_type"
## [9] "maturity_stage" "length_total_cm"
## [11] "Year" "Month"coutlength <-ohbio4 %>%
drop_na(length_total_cm) %>%
dplyr::group_by(Year, Month, asd_code) %>%
dplyr::summarise(avg=mean(length_total_cm))
#kableExtra::kable(coutlength, format = "html")Some features of data set length structures
table(ohbio4$gear_type) %>%
kbl() %>%
kable_styling()| Var1 | Freq |
|---|---|
| Beam Trawls, Midwater | 387680 |
| Gear Not Known, Or Not Elsewhere Included | 43 |
| Midwater Trawls Nei | 14346 |
| Otter Trawls Nei | 3 |
| Otter Trawls, Bottom | 380 |
| Otter Trawls, Bottom, Towed From Stern | 2576 |
| Otter Trawls, Midwater | 375337 |
| Otter Trawls, Midwater, Towed From Side | 1 |
| Otter Trawls, Midwater, Towed From Stern | 1663407 |
table(ohbio2$trawl_technique)%>%
kbl() %>%
kable_styling()| Var1 | Freq |
|---|---|
| C | 436932 |
| T | 2006841 |
Plots
Filter data regarding to previous glances. Follow with a quick glimse to all 48 subarea length composition from monitoring fisheries.
jz <- ggplot(ohbio4 %>%
filter(Year>2000),
aes(x=length_total_cm,
y = as.factor(Year),
fill=asd_code))+
#geom_joy(alpha=0.9) +
geom_density_ridges(stat = "binline", bins = 50,
scale = 1.8,
draw_baseline = FALSE,
alpha=0.9)+
facet_wrap(.~asd_code, ncol=3) +
geom_vline(xintercept = 3.6, color = "red")+
scale_x_continuous(breaks = seq(from = 1, to = 10,
by = 1))+
scale_y_discrete(breaks = seq(from = 2000,
to = 2020, by = 1))+
scale_fill_viridis_d(name="SubArea",
option="F")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90,
hjust = 1))+
xlim(0,10)+
xlab("Longitud (cm.)")+
ylab("")
jz
trwal <- ggplot(ohbio4 %>%
filter(Year>2000),
aes(x=length_total_cm,
y = as.factor(Year),
fill= trawl_technique))+
#geom_joy(alpha=0.9) +
geom_density_ridges(stat = "binline", bins = 50,
scale = 1.8,
draw_baseline = FALSE,
alpha=0.5)+
facet_wrap(.~asd_code, ncol=3) +
geom_vline(xintercept = 3.6, color = "red")+
scale_x_continuous(breaks = seq(from = 1, to = 10,
by = 1))+
scale_y_discrete(breaks = seq(from = 2000,
to = 2020, by = 1))+
scale_fill_viridis_d(name="Trawl Technique",
option="F")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90,
hjust = 1))+
xlim(0,10)+
xlab("Longitud (cm.)")+
ylab("")
trwal
Same plot by Sub Area and month trought year. With this kind of plot, we
can see the recruit power interannually.
ymasd <- ggplot(ohbio4 %>%
filter(Year>2000),
aes(x=length_total_cm,
y = as.factor(Month),
fill=asd_code))+
#geom_joy(alpha=0.9) +
geom_density_ridges(stat = "binline", bins = 50,
scale = 1.5,
draw_baseline = FALSE,
alpha=0.5)+
facet_wrap(.~season_ccamlr, ncol=5) +
geom_vline(xintercept = 3.6, color = "red")+
scale_x_continuous(breaks = seq(from = 1, to = 10,
by = 1))+
scale_y_discrete(breaks = seq(from = 1,
to = 12, by = 1))+
scale_fill_viridis_d(name="SubArea",
option="F")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90,
hjust = 1))+
xlim(0,10)+
xlab("Longitud (cm.)")+
ylab("")
ymasd
yma481 <- ggplot(ohbio4 %>%
filter(Year>2000,
asd_code==481),
aes(x=length_total_cm,
y = as.factor(Month),
fill=asd_code))+
#geom_joy(alpha=0.9) +
geom_density_ridges(stat = "binline", bins = 50,
scale = 1.5,
draw_baseline = FALSE,
alpha=0.5)+
facet_wrap(.~Year, ncol=5) +
geom_vline(xintercept = 3.6, color = "red")+
scale_x_continuous(breaks = seq(from = 1, to = 10,
by = 1))+
scale_y_discrete(breaks = seq(from = 1,
to = 12, by = 1))+
scale_fill_viridis_d(name="SubArea",
option="F")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90,
hjust = 1))+
xlim(0,10)+
xlab("Longitud (cm.)")+
ylab("")
yma481
Mean Length
Recruit estimate (% Low Recruit length) (Perry, 2020)
Un diaggrama de caja para el porcentage de individuos bajo talla (3.6 mmm.)
box <- ggplot(ohbio4) +
geom_boxplot(aes(length_total_cm, group=season_ccamlr),
alpha=0.3,
fill=3)+
xlim(0, 7)+
coord_flip()+
facet_wrap(~asd_code, ncol=1)+
theme_bw()+
labs(x="% Below Recruit length (3.6 mm.)",
x="")
hi <- ggplot(ohbio4)+
geom_histogram(aes(y=length_total_cm),
fill=3,
alpha=0.3,
color="black")+
coord_flip()+
theme_bw()+
ylim(0,10)
box | hi
Ahora saco tallas medias por distintas variables
pmea <- ggplot(coutlength %>%
filter(Year>2000),
aes(Year,avg))+
geom_point(shape=21, fill=coutlength$asd_code,
show.legend = T) +
stat_smooth(method= "glm", colour='#253494')+
scale_size(range = c(-4,8)) +
theme_bw()+
facet_wrap(.~asd_code)+
scale_x_continuous(breaks = seq(from = 2000, to = 2020, by = 2))+
#scale_y_discrete(breaks = seq(from = 1, to = 13, by = 1))+
theme(axis.text.x = element_text(angle = 90, hjust = 2))+
guides(fill = guide_legend(reverse=F))+
scale_fill_viridis_c(option="E")+
ylim(3,6)+
ylab("") +
xlab("") +
ggtitle("Lenght Mean Krill fishery")
pmea
by month and year only 48.1
pmea481 <- ggplot(coutlength %>%
filter(Year>2000,
asd_code==481),
aes(Month, avg))+
geom_point(shape=21,
show.legend = T) +
stat_smooth(method= "glm", colour='#253494')+
scale_x_continuous(breaks = seq(from = 1, to = 12, by = 2))+
theme_bw()+
facet_wrap(.~Year) +
ggtitle("Lenght Mean Krill fishery 48.1")+
ylim(0, 6)
pmea481
GLM
model1 <- glm(avg ~ Year, data = coutlength)
r2(model1)## R2: 0.000check_model(model1)
model2 <- glm(avg ~ Year +
Month, data = coutlength)
r2(model2)## R2: 0.009check_model(model2)
model3 <- glm(avg ~ Year +
Month +
asd_code, data = coutlength)
r2(model3)## R2: 0.013check_model(model3)
Compara modelos
plot(compare_performance(model1,
model2,
model3,
rank = TRUE))
Probar datos con RaadTools library (Sumner, 2022)
References
Sea ice data is from Giovanni NASA browser (https://giovanni.gsfc.nasa.gov/) ref: Global Modeling and Assimilation Office (GMAO) (2015),
MERRA-2 tavgM_2d_flx_Nx: 2d,Monthly mean,Time-Averaged,Single-Level, Assimilation,Surface Flux Diagnostics V5.12.4, Greenbelt, MD, USA, Goddard Earth Sciences Data and Information Services Center (GES DISC), Accessed: [06 July 2022], 10.5067/0JRLVL8YV2Y4